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The static and dynamical properties of heavy water have been studied at ambient conditions 
with extensive Car-Parrinello molecular-dynamics simulations in the canonical ensemble, with tem- 
peratures ranging between 325 K and 400 K. Density-functional theory, paired with a modern 
exchange-correlation functional (PBE), provides an excellent agreement for the structural proper- 
ties and binding energy of the water monomer and dimer. On the other hand, the structural and 
dynamical properties of the bulk liquid show a clear enhancement of the local structure compared 
to experimental results; a distinctive transition to liquid-like diffusion occurs in the simulations only 
at the elevated temperature of 400 K. Extensive runs of up to 50 picoseconds are needed to obtain 
well-converged thermal averages; the use of ultrasoft or norm-conserving pseudopotentials and the 
larger plane-wave sets associated with the latter choice had, as expected, only negligible effects on 
the final result. Finite-size effects in the liquid state are found to be mostly negligible for systems 
as small as 32 molecules per unit cell. 



I. INTRODUCTION 

Water, due to its abundance on the planet and its role 
in many of the organic and inorganic chemical processes, 
has been studied extensively and for decades both at the 
theoretical and at the experimental leveli - — . The pecu- 
liar interplay of hydrogen bonding, glassy behavior, and 
of quantum-mechanical effects on the dynamics of the 
atomic nuclei make computer simulations challenging, 
and a great effort has been expended to build a com- 
prehensive and consistent microscopic picture, and a link 
with observed macroscopic properties' - ^!. Additionally, 
it is only recently that close agreement for fundamen- 
tal structural information such as the radial distribution 
function has been obtained between different experimen- 
tal techniques, such as X-raj*^ and neutron diffraction 3 
measurements. 

Computational studies based on molecular dynamics 
simulations have also a rich history in the field. Simu- 
lations using force-fields models' - — have been successful 
at reproducing many structural and dynamical proper- 
ties of liquid water. However, empirical models rely on 
parameters which are determined by fits to known exper- 
imental data, or occasionally to ab-initio results. Their 
transferability to different environments, or the ability 
to reproduce faithfully the microscopic characteristics of 
hydrogen bonding, are often in question. Due to develop- 
ment of novel techniques"^ and the ever-increasing im- 
provement in computational power, extensive molecular- 
dynamics simulations from first-principles are now possi- 
ble. The increased accuracy and predictive power of these 
simulations comes at a significant price, and careful con- 
siderations has to be given to the length scales and time 
scales that can be afforded in a first-principle simulation, 
and the trade-offs in statistical errors when compared 



with classical simulations. Numerous ab-initio simula- 
tions on water have appeared^ - — , showing good agree- 
ment with experiments for the structural and dynamical 
data. Recent results have also reported2£*2i that careful 
equilibration for ab-initio water at ambient conditions 
leads to radial distribution functions over-structured in 
comparison with experiments 2 . After equilibration, the 
numerical estimates for the diffusion coefficient become 
at least one order of magnitude smaller than the mea- 
sured ones. 



Prompted by these results and by our own observations 
on the structure of water around iron aqua ions 2 ^, we 
have undertaken an extensive investigation of the static 
and dynamical properties of water, to ascertain its phase 
stability around ambient conditions as predicted by first- 
principles molecular dynamics. Particular care has been 
given to the statistical accuracy of the results, assuring 
that the time scales and length scales of the simulations 
were chosen appropriately for the given conditions. The 
paper is organized as follows: In section II, we detail all 
the technical aspects of our simulations. Section III sur- 
veys the static and vibrational properties of the water 
molecule and the water dimer in vacuum, at the GGA- 
PBE 23 density-functional level. In section IV, we discuss 
our extensive liquid water simulations, performed with 
Car-Parrinello molecular dynamics, in the temperature 
range between 325 K and 400 K. Section V discusses the 
limitations of this approach, and some of the possible 
reasons to explain the remaining discrepancies with ex- 
perimental results. 



2 



II. CAR-PARRINELLO MOLECULAR 
DYNAMICS 

Our first-principles calculations are based on density- 
functional theory, periodic-boundary conditions, plane- 
wave basis sets, and norm-conserving24 or ultrasoft 
pseudopotentials 2 ^ to represent the ion-electron interac- 
tions, as implemented in the public domain codes CP and 
PWSCF in the ^-ESPRESSO package^. 

In Car-Parrinello molecular dynamics, an extended La- 
grangian is introduced to include explicitly the wavefunc- 
tion degrees of freedom, that are evolved "on-the-fly" si- 
multaneously with the ionic degrees of freedom: 

L C p = dr |* i(r) f + 5^ Ml ^~ 

i I 

E KS [{*i},Rj] + 

Y^Ki^j dv%{v)^ j {v)-8 i ^j. (1) 

In Eq. ^ (r) are the occupied Kohn-Sham orbitals, /i 
is a fictitious mass parameter used to control the evolu- 
tion of the electronic degrees of freedom in time, Mj are 
ion masses, Eks is the Kohn-Sham energy and Ay are 
Lagrange multipliers, used to impose the orthonormal- 
ity constraint J = Sij. An extensive review of the 

method can be found in Refill; the subtle technical issues 
arising in the simulations of liquid water have been dis- 
cussed in the recent literatur e 2 ? i 2 ? . We note in passing 
that our simulations' parameters agree with the recom- 
mendations set forth in these latter papers. It should be 
stressed that full convergence of the forces acting on the 
nuclei with respect to the basis set is easily achieved with 
plane waves, and no Pulay forces arise in the dynamics, 
since the basis set is independent from the atomic posi- 
tions. 



A. Technical Details 

The structural and vibrational properties of the wa- 
ter monomer and dimer and the binding energy of 
the dimer have been calculated using density-functional 
theory in the generalized-gradient approximation and 
the total energy pseudopotential method, and density- 
functional perturbation theory^, as implemented in 
PWSCF. We performed separate calculations using ei- 
ther norm-conserving pseudopotentials for both the hy- 
drogen and the oxygen, or ultrasoft ones. These same 
pseudopotentials were also used for the norm-conserving 
or ultrasoft molecular dynamics simulations. In particu- 
lar, the O Troullier-Martins norm-conserving pseudopo- 
tential was generated using the FHI98PP package^ , with 
core radii for the s, p and d components of 1.25 a.u., 1.25 
a.u., and 1.4 a.u. respectively. The Troullier-Martins 
hydrogen pseudopotential was generated using the Atom 
code^i with a core radius for the s component of 0.8 



TABLE I: Structural properties of the water monomer and 
dimer and binding energy of the dimer, as obtained in DFT- 
PBE using ultrasoft or norm-conserving pseudopotentials, 
and compared to available experimental and theoretical re- 
sults. 





PBE US 

(This 

work) 


PBE NC 

(This 

work) 


PBE NC 
(Re^ 




BLYP^i 


t-HOH 


104.6 U 


104.2 U 


104.2" 


104.5" 


104.4" 


doif(A) 


0.98 


0.97 


0.97 


0.96 


0.97 


Loho 


173° 


172° 


174° 


174° 


173° 


doo(A) 


2.89 


2.88 


2.90 


2.98 


2.95 


^dimer 


-23.2 


-23.8 


-21.4 


-22.8 


-18 


(kJ/mol) 













a.u. The ultrasoft pseudopotentials were taken from the 
standard PWSCF distribution^. The Kohn-Sham or- 
bitals and charge density have been expanded in plane 
waves up to a kinetic energy cutoff of 25 and 200 Ry (re- 
spectively) for the ultrasoft case, and of 80 Ry and 320 
Ry for the norm-conserving case. A cubic supercell of 
side 30 a.u. was used; interaction with periodic images 
is negligible 3 ^ with this unit cell size. 



III. WATER MONOMER AND DIMER: 
STRUCTURAL AND VIBRATIONAL 
PROPERTIES 

The equilibrium structures and energetics are summa- 
rized in Tabled We have included published results^ us- 
ing the BLYP functional for comparison. Both ultrasoft 
and norm-conserving PBE density functionals show very 
good agreement with experimental values. In particular, 
the PBE results have a dimer binding energy in closer 
agreement to the experiments than BLYP; the binding 
energy in this latter case is too weak by 4 kJ/mol, and 
exhibits a longer 0-0 distance. 

Table ITU and lllll show respectively the vibrational fre- 
quencies of the water monomer and dimer in vacuum. In 
this calculation, a hydrogen mass of 1 a.m.u. was used 
(as opposed to the 2 a.m.u. mass used for the dynami- 
cal simulations of heavy water). The calculations of the 
vibrational frequencies with PWSCF using a cubic cell 
of size (30 a.u.) 3 . To achieve a convergence of a few 
cm -1 in the frequencies, cutoffs of 35 Ryd and 420 Ryd 
were used for wavefunctions and charge densities, respec- 
tively, in the ultrasoft case and, in the norm-conserving 
case, 100 Ryd and 400 Ryd were used. As shown in Ta- 
ble HJ and Unlthe PBE functional gives intramolecular 
stretching modes that are in general blue-shifted com- 
pared to BLYP and experimental results. In the calcula- 
tions for the dimer, the libration modes are also higher 
for the PBE functionals then those given by experiments 
and BLYP. We note in passing that the errors on these 
frequencies (especially the low energy ones) are slightly 
larger than usually expected from DFT). We will return 
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TABLE II: Vibrational frequencies of water monomer: v\ , 1/2 
and V3 are the symmetric stretching, bending and asymmetric 
modes, respectively. 





PBE (US) 


PBE (NC) 


Expt^S 


BLYP — 




3781 


3704 


3657 


3567 




1573 


1599 


1595 


1585 


^ 3 (cm _1 ) 


3908 


3816 


3756 


3663 



TABLE III: Vibrational frequencies of water dimer: v\, V2 
and axe the symmetric stretching, bending and asymmet- 
ric stretching modes, respectively. Proton acceptor and donor 
molecules are denoted as (A) and (D). v(Hb) are the two libra- 
tion modes between molecules and v(0 — O) is the hydrogen- 
bond stretching mode. 





PBE (US) 


PBE (NC) 


Expt 34 > 39 


BLYP 33 




3778 


3695 


3622 


3577 


^(Ajfcm" 1 ) 


1570 


1596 


1600 


1593 


v.i{A)(cm~ 1 ) 


3901 


3804 


3714 


3675 


i/i(fl)(cm _1 ) 


3601 


3532 


3548 


3446 


u 2 (D)(cm' 1 ) 


1593 


1616 


1618 


1616 


i/s(.D)(cro -1 ) 


3871 


3781 


3698 


3647 


i/tffljfcm" 1 ) 


666 


644 


520 


600 


!/(#&) (cm -1 ) 


379 


378 


320 


333 


u(0 - 0)(cmT x ) 


202 


196 


243 


214 



to this point in later section. 

IV. LIQUID WATER SIMULATIONS 
A. Liquid water simulation at 325 K 

1. Simulation Details 

In this first simulation, we used a body-centered-cubic 
supercell with 32 heavy water molecules, periodic bound- 
ary conditions, and the volume corresponding to the 
experimental^ density of 1.0957 g/cm 3 at 325 K. A 
body-centered-cubic supercell strikes the optimal bal- 
ance, for a given volume, in the distance between a 
molecule and its periodic neighbors, and the number 
of these periodic neighbors. Ultrasoft pseudopotentials 
were first used, as detailed in the previous section, with 
plane-wave kinetic energy cutoffs of 25 Ry (wavefunc- 
tions) and 200 Ry (charge densities). The deuterium 
mass was used in place of hydrogen to allow for a larger 
time step of integration. It should be noted that for clas- 
sical ions this choice does not affect thermodynamic prop- 
erties such as the melting temperature (the momentum 
integrals for the kinetic energy factor out in the Boltz- 
mann averages). Of course, dynamical properties such as 
the diffusion coefficient will be affected by our choice of 
heavier ions. Extensive experimental data for deuterated 
(heavy) water are in any case widely available. 

The wavefunction fictitious mass (/z) is chosen to be 



700 a.u.; this results in a factor of ~14 between the aver- 
age kinetic energy of the ions and that of the electrons. 
A time step of 10 a.u. was used to integrate the electron 
and ionic equations of motions. This combined choice of 
parameters allows for roughly 25 ps of simulation time 
without a significant drift in the kinetic energy of wave- 
functions and the constant of motion for the Lagrangian 
(JTJ. Our choice of fictitious mass is consistent with the 
ratio fi/M < | for heavy water molecules suggested by 
Grossman et al— , and assures that the physical proper- 
ties are not influenced by the electronic degrees of free- 
dom. Our initial configuration was obtained from a com- 
paratively short 1.2 ps simulation at twice the value of 
the target temperature (650 K instead of 325 K); a restart 
with zero initial velocities was then performed at 325 K, 
with the temperature controlled by a single Nose-Hoover 
thermostat on the ions (no electronic thermostat is ap- 
plied in any of the simulations). 



2. Results 

The thermostat stabilizes quickly (~1 ps) the temper- 
ature around the target value of 325 K, with the usual 
fluctuations due to small size of the system. However, the 
system is still far from equilibrium; this can be clearly 
observed by looking at the time evolution of the radial 
distribution function (RDF) and mean square displace- 
ments (MSDs). We plot in Fig.[T]the MSDs of the oxygen 
atoms as a function of time. For the first 10 ps, the water 
molecules diffuse with a velocity comparable with exper- 
imental data; after 10 ps, a sharp drop in diffusivity is 
observed, accompanied by a distinctive sharpening of the 
features in the oxygen-oxygen radial distribution function 
(see Fig.|2J|. The radial distribution functions was calcu- 
lated from the infinite bulk system by repeating the unit 
cell in all directions. This means that molecules up to 
5.5 A in a 32-molecule cell are inequivalent. As we move 
beyond this cutoff distance, the radial distribution func- 
tions will include both molecules that are inequivalent, 
and some that are equivalent. The statistical accuracy is 
going to be gradually, but slowly, affected. In this graph, 
and most of the latter figures, we plot radial distribution 
functions only up to 5.5 A(but see e.g. Fig. 10 for a 
discussion of the finite-size effects). 

The potential energy for the dynamics (Fig. also 
drifts downward in these first 10 ps, and stabilizes after- 
ward. We calculated the self-diffusion coefficient (D se if) 
from the Einstein relation (in 3 dimensions): 

6^/ = ^™ ^(|r 4 (t)-r,(0)| 2 ). (2) 

The structural and dynamical properties before and after 
this 10 ps mark are summarized in Table llVl experimen- 
tal values at 298 K are included for comparison. All 
these observation conjure to a picture in which the sys- 
tem takes at least 10 ps to reach a reasonably thermalized 
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FIG. 1: Mean square displacement and potential energy as a 
function of time for our first 32- water molecules simulations at 
325 K. A: Diffusive region in the first stage of the simulation. 
B: After about 10 ps diffusivity drops abruptly. 



FIG. 2: O-O radial distribution function before (A) and after 
(B) the equilibrium is attained. 



TABLE IV: Structural and dynamical parameters before and 
after the fO ps mark, compared with the experimental results 
at 298 K. 







D (cm 2 /s) 


Before 


2.82 


1.5 x l(T b 


After 


3.21 


0.14 x 10 -5 


Expfeffi^ 


2.75 


2.0 x 10~ 5 



state, in a process somewhat reminiscent of a glass tran- 
sition. Although the time needed for equilibration will 
be dependent on the initial conditions, these preliminary 
result suggests that simulation times in the order of ten 
of picoseconds might be need to calculate well-converged 
thermodynamic observables. 

Once the initial thcrmalization trajectory was dis- 
carded from our averages, we obtained a self-diffusion 
coefficient one order of magnitude smaller than the ex- 
perimental value measured at room temperature. This 
result, combined with the clear over-structuring of the 
oxygen-oxygen radial distribution function goo( r )i indi- 
cates that our system has reached a "frozen" equilibrium 
state very different from what expected for liquid water 
(strictly speaking, a system with a finite and small num- 
ber of incquivalcnt atoms or molecules will never undergo 
a phase transition). 

These considerations indicates the need to assess accu- 
rately the phase stability of liquid water as obtained from 
first-principle molecular dynamics. At the same time, 
they point out the requirement of long simulation times, 
and a careful analysis of the technical details of the sim- 
ulations. 



0.4 




ps 



FIG. 3: Mean square displacements as a function of time of 
our initial equilibration runs for 32 water molecules at 325, 
350, 375 and 400 K. 



B. Extensive water simulations in the region 
between 325 K and 400 K 

1. Simulation details 

In order to find out the temperature at which our sys- 
tem would move from a glassy to a liquid-like state, we 
decided to perform a series of extensive simulations at 
increasing temperature, from 325 K to 400 K. We first 
performed w 25 ps simulations at 325 K, 350 K, 375 K 
and 400 K, using /x=700 a.u. and St=10 a.u. . We used 
at every temperature the experimental densities^, with 
the caveat that the 400 K value was obtained by extrapo- 
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FIG. 4: Potential energy and constant of motion in a produc- FIG. 6: 0-0 radial distribution functions calculated from the 
tion run at 400 K. first and the next 12 ps of the simulation at 400 K. 
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FIG. 5: Kinetic energy of the ions and the electrons in a 
production run at 400 K. 



lation. 25 ps is approximately the maximum time allowed 
for a simulation with these parameters before the drift in 
the kinetic energy of the wavefunctions becomes appar- 
ent. We thus used these simulations as efficient "ther- 
malization" runs, to be followed by production runs that 
will be described below. It is interesting to monitor dur- 
ing these thermalization runs the evolution of the MSDs; 
these are shown for all four temperatures in Fig. [3J An 
abrupt drop in diffusivity is observed for all cases but 
one (predictably, the one at the highest temperature of 
400 K). The onset of this drop in diffusivity varies, but 
broadly speaking is again of the order of 10 ps. 

With these trajectories, now well thermalized in con- 
figuration space, we started our four production runs at 



325 K, 350 K, 375 K and 400 K, and each of them starting 
from the last ionic configurations of the previous simu- 
lation at the corresponding temperature, but with zero 
ionic velocities. These four runs, lasting between 20 ps 
and 37 ps, were performed using /i=450 a.u. and 5t=7 
a.u. This choice of mass and timestep allows for an ex- 
cellent conservation of the constant of motion, and negli- 
gible drift in the fictitious kinetic energy of the electrons, 
for the simulation times considered. The ratio between 
the kinetic energy of the ions and that of the electrons 
was « 22 for the whole production time. Although a 
small enough fictitious mass decouples the electronic and 
ionic degrees of freedom, Tangney et. al4i pointed out 
that there is a fictitious mass dependent error that is 
not averaged in the time scale of ionic motions. Schwc- 
gler et. ali^S studied this effect comparing closely Car- 
Parrinello and Born-Oppenheimer simulations finding a 
larger self-diffusion coefficient in the Car-Parrinello sim- 
ulation. However, the structural and thermodynamical 
properties were not affected. 

We show in Figs. 01 and [5] the case of the 400 K simu- 
lation; we stress that no periodic quenching of the elec- 
trons was needed, and the simulations were single un- 
interrupted runs. Since the initial configurations were 
already at equilibrium at their respective temperature, 
and thermalization in momentum space is fast, we found 
that "production time" can start early in the simulations. 
We discarded from each trajectory the initial 1.2 ps that 
were needed to allow the ions to reach their target ki- 
netic energy. As a measure of the good thermalization 
reached in the simulations, we show in Fig. [S] the 0-0 
radial distribution function obtained from the first 12 ps 
of our 400 K trajectory, and the following 12 ps. 

To rule out any spurious effect in our simulations com- 
ing from the use of pseudopotentials, or an extended La- 
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TABLE V: Details of the production runs. 



Pseudo- 


TfKl 


Density 




St 


Production 


potentials 




(g/cm 3 ) 


(a.u.) 


(a.u.) 


time (ps) 


US 


325 


1.0957 


450 


7 


37.6 


NC 


325 


1.0957 


300 


5 


22.1 


US 


350 


1.0815 


450 


7 


22.9 


us 


375 


1.0635 


450 


7 


21.1 


us 


400 


1.0554 


450 


7 


32.5 


NC 


400 


1.0554 


300 


5 


20.2 



grangian, we also performed two simulations using norm- 
conserving pseudopotentials (as described in Section II). 
These require larger plane-wave basis sets (80 Ry for the 
wavefunctions and 320 Ry charge densities, correspond- 
ing to 40000 plane waves vs 7000 for the ultrasoft case), 
and discrepancies, if any, with the ultrasoft calculation 
will provide an approximate estimate of the effects of the 
pseudopotential approximation and of the dynamics of 
the fictitious degrees of freedom. For the runs involving 
these norm-conserving pseudopotentials, we used /i =300 
a.u. and 5t—5 a.u^. These parameters results in a factor 
of ~ 13 between the kinetic energy of the ions and that of 
the electrons. Details of all these simulations are shown 
in Table E| 



2. Results 

The oxygen-oxygen and oxygen-deuterium radial dis- 
tribution functions for the different conditions consid- 
ered in this work are shown in Fig. [7| and Fig. [S] re- 
spectively. We find that at temperatures of 375 K and 
below both goo and goH show considerably more struc- 
ture than found experimentally at 300 K. For this range 
of temperatures the height of the first peak of goo (r) is 
roughly between 3.2 and 3.4, and significantly larger than 
the experimental value of 2.75 (also measured at 300 K). 
However, when the temperature in our simulations is in- 
creased to 400 K, a distinct drop of the first peak to ~ 
2.5 is observed, the radial distribution function goo{f) 
and goHir) show a sharp change in their structure, and 
the water molecules start diffusing much faster, as re- 
flected in the MSDs curves for the oxygen atoms shown 
in Fig. To provide cleaner statistics, the MSDs curves 
shown have been calculated as an average over individ- 
ual MSDs curves, each obtained from our trajectory by 
shifting - for each individual MSDs curve - the start- 
ing configuration by 0.017 ps (in this way, a 17 ps tra- 
jectory would provide 1000 progressively shorter MSDs 
curves that are then averaged). The self-diffusion coeffi- 
cient D se i f is calculated from the slope of the respective 
MSDs curve in the range of 1 to 20 ps using Einstein's 
relation ©. Negligible differences are observed between 
simulations performed with ultrasoft or norm-conserving 
pseudopotentials, ruling out any role of the pseudopoten- 
tial details in this observed result. 



TABLE VI: Summary of structural and dynamical properties 
of water. D se if is self-diffusion coefficient. g max is first peak 
height and R[j moI ] is location of first peak. The last column 
is the number of Hydrogen-bond per molecule. 





Dseif 

10" 5 cm 2 /s 


Qmax 




no. of H-bonds 
per molecule 


325 K (US) 


0.07 


3.38 


2.67 


3.86 


325 K (NC) 


0.16 


3.25 


2.73 


3.79 


350 K (US) 


0.25 


3.25 


2.73 


3.77 


375 K (US) 


0.26 


3.10 


2.71 


3.72 


400 K (US) 


2.03 


2.50 


2.73 


3.45 


400 K (NC) 


1.66 


2.55 


2.75 


3.37 


Expl42d& at 300 K 


1.80 


2.75 


2.80 


3.58 



The structural and dynamical results are summarized 
in Table IVT1 As seen in this table, there is an eight-fold 
increase in D se if when increasing the temperature from 
375 K to 400 K. Price et. al. 43 reported the experimen- 
tal self-diffusion coefficient of supercooled heavy water at 
different of temperatures. At 276.4 K, which is just be- 
low the freezing temperature of heavy water (277.0 K), 
the experimental value for D se if is 0.902 x 10~ 5 cm 2 /s. 
D se if for our simulations at 325 K, 350 K and 375 K 
is 0.16, 0.25 and 0.26 xl0 _5 cm 2 /s, respectively. These 
numbers are significantly smaller then the experimental 
value below the freezing point. On the other hand, D se [f 
at 400 K in our simulations is comparable to the exper- 
imental value at 300 K. These observations suggest that 
the theoretical freezing point for water at the DFT-PBE 
level is between 375 K and 400 K, and water below 375 
K is in a glassy/supercooled state. 

The hydrogen-bond structure can be studied calculat- 
ing the number of hydrogen bonds per molecule: we iden- 
tify a hydrogen bond is identified when two oxygen atoms 
are closer than 3.5 A and the Loho angle is greater than 
140° (consistently with Re£^i, and at slight variance with 
ReS% The results are shown in the last column of Ta- 
ble E] Between the temperatures of 325 K and 375 K 
there are only small changes in the number of hydrogen 
bonds in the system, going from 3.86 per molecule to 3.72 
per molecule. An abrupt decrease to 3.37-3.45 bonds oc- 
curs when the temperature increases from 375 K to 400 
K. The experimental value of 3.58 at 300 K lies between 
our values at 375 K and 400 K, further suggesting that 
the freezing point in our simulations is between 375 K 
and 400 K. 

In summary, we found clear liquid-like signatures in the 
structure and dynamics of (heavy) water, as described 
by DFT-PBE and Car-Parrinello MD, for temperatures 
reaching at least 400 K. At temperatures of 375 K and 
lower water is found to be in a glassy state, more struc- 
tured and with much lower diffusivity. This discrepancy 
of more than 100 K between experimental and theoretical 
results is obviously relevant, given the enormous impor- 
tance of water in the description of systems ranging from 
electrochemistry to biology, and it is investigated further 
in the next section. 
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Radius (A) 



FIG. 7: 0-0 radial distribution functions calculated from the 
production runs at 325 K, 350 K, 375 K and 400 K for ul- 
trasoft and norm-conserving pseudopotentials. Experimental 
result is taken from ReS». 



FIG. 9: Mean square displacements calculated from simula- 
tions at 325 K, 350 K, 375 K and 400 K for ultrasoft and 
norm-conserving pseudopotentials. 



A. Finite-size effects 
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FIG. 8: O-D radial distribution functions calculated from 
simulations at 325 K, 350 K, 375 K and 400 K for ultrasoft 
and norm-conserving pseudopotentials. Experimental result 
is taken from Rei-. 



OVERESTIMATION OF THE FREEZING 
TEMPERATURE OF WATER 



There are several possible reasons for the overestima- 
tion of the freezing temperature of water, and several of 
them could play a significant role. We discuss here some 
of the possibilities. 



Since our simulation cell contains only 32 water 
molecules, finite-size effects could obviously play a role 
even if periodic-boundary conditions are used. The in- 
teractions of water molecules with their periodic images 
could be considerable due to the long-range hydrogen- 
bond network. On the other hand, when zero correla- 
tions are found between a molecule and its periodic image 
we can safely assume that the unit cell is for all practi- 
cal purposes large enough, and every molecule feels the 
same environment that it would have in an infinite sys- 
tem. In our case, the distance between a molecule and its 
eight periodic images is ~ 11 A and at this distance all 
radial distribution functions look very flat and unstruc- 
tured. In any case, to study the finite-size effects we car- 
ried out another extensive simulation (40 ps total, with 
15 ps of production time following 25 ps of thermaliza- 
tion) for a system composed of 64 heavy water molecules 
at 400 K. We used the same parameters for this simu- 
lation as in the 32-molcculc, 400 K ultrasoft simulation. 
The oxygen-oxygen radial distribution function goo (r) is 
shown in Fig. ^| for both the 32- and 64-molecule sys- 
tems. As mentioned previously, the radial distribution 
functions was calculated by repeating the unit cell in all 
directions. The molecules up to 5.5 A in the 32-molecule 
cell, and 6.9 A in the 64-molecule cell are inequivalent. 
We indicate in the graphs with two arrows, the radii of 
the spheres completely inscribed by in BCC simulation 
cells. 

Differences between 32- and 64-molecule systems are 
negligible, and within the variance for simulations of the 
order of 10-20 ps (as estimated from uncorrelated clas- 
sical simulation data™); the 64- water simulation shows 
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FIG. 10: O-O radial distribution function for a Car-Parrinello 
simulations with 32 or 64 molecules. The arrows indicate the 
radius of the sphere that is completely inscribed by the BCC 
simulation cells with 32 or 64 molecules. 



FIG. 12: Power spectrum of deuterium atoms calculated from 
the velocity-velocity correlation function. For comparison, 
we also show the power spectrum as obtained with a larger 
fictitious mass of 700 a.u (dotted line) instead of 450 a.u. 
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FIG. 11: O-O radial distribution function for a classical (SPC) 
simulation with 64 or 1000 water molecules. 



a marginally more structured g(r) where the first peak 
height is at about 2.6 (compared to 2.5 for the 32-water 
case) . Larger ab-initio simulations would be too demand- 
ing; for this reason, we performed two classical simu- 
lations at 300 K using the SPC force field^ for water, 
and comparing the case of 64 and 1000 water molecules 
(each simulation lasting 1000 ps). The goo{i) calculated 
from these two runs are shown in Fig. ^2 and again we 
do not find any significant differences between these two 
curves. These results help ruling out finite-size effects 
as the major cause of the discrepancy observed with the 
experimental numbers. 



B. Exchange-correlation functional effects 



While density- functional theory is in principle exact, 
any practical application requires an approximated guess 
to the true exchange-correlation functional. In this work, 
we have used the GGA-PBE approach 23 . As it was ob- 
served in Sec. Ill, the structural properties for the water 
molecule and dimer are in excellent agreement with ex- 
periments, as is the binding energy for the dimer. On 
the other hand, the vibrational properties show larger 
discrepancies with experiments than usually expected, in 
particular for some of the libration modes in the dimer. 
This result certainly points to the need for improved func- 
tionals to describe hydrogen bonding. The dependence 
of the melting point on the exchange-correlation func- 
tional chosen is more subtle; below 400 K, PBE water dis- 
plays solid-like oxygen-oxygen radial distribution func- 
tions that are only slightly affected by the temperature, 
and that are similar to those obtained with a fairly differ- 
ent functional such as BLYR5 -. The similarity between 
these radial distribution functions is just a reminder that 
the structural property and the geometry of the inter- 
molecular bonds are well described by different function- 
als; once water is "frozen" , all radial distribution func- 
tions will look similar. The temperature at which this 
transition takes place could be affected by the use of dif- 
ferent functionalist and the magnitude of the contribu- 
tion of any one of them to the melting point temperature 
is still open to investigation. 
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C. Quantum effects 

In first-principles Car-Parrincllo or Born-Oppcnhcimcr 
molecular dynamics simulations the ions are most of- 
ten treated as classical particles, which is a good ap- 
proximation for heavy ions (path-integral simulations 
can describe the quantum nature of the ionsa^ - — , but 
their computational costs, when paired with a first- 
principles DFT descriptions of the electrons, preclude 
at this moment simulations with the statistical accuracy 
needed). However, for light ions like hydrogen or deu- 
terium, the effects of a proper quantum statistics can be 
very significant; tunneling of the nuclei can also affect the 
dynamics^^. In the case of water, all the intramolec- 
ular vibrational modes and some of the intermolecular 
modes are much higher in energy compared to room tem- 
perature. We show in Fig. ll2l the power spectrum for the 
deuterium atoms as calculated from the velocity-velocity 
correlation function of heavy water molecules for our sim- 
ulations at 400 K (ultrasoft, 32 molecules). The distinc- 
tive peaks of the intramolecular stretching and bending 
modes are centered around 2400 cm' 1 and 1200 cm -1 
respectively, much higher than the room temperature of 
ksTroom ~ 200 cm -1 . The peak at 500 cm" 1 corresponds 
to the intermolecular vibrational modes, also larger than 
kBT room . When ions are treated as classical particles, 
as in our ab-initio molecular dynamics simulations, all 
vibrational modes obey Boltzmann statistics. In reality, 
modes with frequency higher than fcsT room are frozen 
in their zero-point motion state, and their exchange of 
energy with the lower-frequency modes ("the environ- 
ment") is suppressed - in other words their contribution 
to the specific heat is zero, in full analogy with the low- 
temperature discrepancies from the Dulong-Petit law in 
solids. This effect could significantly affect the dynamics 
of water-water interactions, and it has long been argued 
that treating each water molecule as a rigid body could 
actually provide a closer match with experimental con- 
ditions. In fact, recent ab-initio simulations** in which 
the water molecules are constrained to maintain their 
equilibrium intramolecular bond lengths and bond angle 



result in a more diffusive and less structured description 
of water, that remains liquid at a temperature of 326 K. 
Path-integral simulation** for water described with clas- 
sical potentials did find as well a significant difference 
due to the quantum effects (i.e. the freezing of the high- 
energy vibrational excitation), of the order of 50K. 



VI. CONCLUSION 

We performed extensive first-principles molecular dy- 
namics simulations of heavy water at the DFT-PBE level. 
Equilibration times are found to be comparatively long, 
and easily in excess of ks 10 ps at ambient temperature. 
Good statistics was obtained with simulations that in- 
cluded at least 25 ps of thermalization time, followed by 
20 to 40 ps of production time. At ambient tempera- 
ture, water is found to be over-structured compared to 
experiment, and with a diffusivity that is one order of 
magnitude smaller than expected. An abrupt change in 
the structure and dynamics is observed when the the tem- 
perature is raised from 375 K to 400 K; even at this high 
temperature, where a liquid-like state is reached, water 
shows more structure and less diffusivity than found ex- 
perimentally at room temperature. These results are 
broadly independent of some of the possible errors or 
inaccuracies involved in first-principles simulations, in- 
cluding in this case finite-size effects, insufficient ther- 
malization or simulation times, or poorly designed pseu- 
dopotentials. Our simulations suggest that the freezing 
point is around 400 K - this discrepancy is at variance 
with the good agreement for the structure and energetics 
of the water monomer and dimer with the experimen- 
tal values, and could originate in the neglect of quantum 
statistics for the many high-frequency vibrational modes 
in this system. 

This work was performed under the support of the 
Croucher Foundation and Muri Grant DAAD 19-03-1- 
0169. Calculations in this work have been done using 
the ^-ESPRESSO package** 



1 Water: A comprehensive Treatise, edited by F. Franks 
(Plenum, New York, 1972) Vol 1. 

2 J. M. Sorenson, G. Hura, R. M. Glaeser, and T. Head- 
Gordon, J. Chem. Phys. 113, 9149 (2000). 

3 A. K. Soper, Chem. Phys. 258, 121 (2000). 

4 A. Rahman and F. H. Stillinger, J. Chem. Phys. 55, 3336 
(1971). 

5 W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. 
Impey, and M. L. Klein J. Chem. Phys. 79, 926 (1983). 

6 K. Toukan and A. Rahman, Phys. Rev. B 31, 2643 (1985). 

7 S. W. Rick, S. J. Stuart, and B. J. Berne, J. Chem. Phys. 
101, 6141 (1994). 

8 G. Lamoureux, A. D. MacKerell, Jr., and B. Roux, J. 
Chem. Phys. 119, 5185 (2003). 



S. Izvekov, M. Parrinello, C. J. Burnham and A. Voth, J. 
Chem. Phys. 120, 10896 (2004). 

10 K. Laasonen, M. Sprik, and M. Parrinello, J. Chem. Phys. 
99, 9080 (1993). 

11 M. Bernasconi, P. L. Silvestrelli, and M. Parrinello, Phys. 
Rev. Lett. 81, 1235 (1998). 

12 P. L. Silvestrelli and M. Parrinello, Phys. Rev. Lett. 82, 
3308 (1999). 

13 P. L. Silvestrelli and M. Parrinello, J. Chem. Phys. Ill, 
3572 (1999). 

14 E. Schwegler, G. Galli, and F. Gygi, Phys. Rev. Lett. 84, 
2429 (2000). 

15 M. Boero, K. Terakura, T. Ikeshoji, C. C. Liew, and M. 
Parrinello, J. Chem. Phys. 115, 2219 (2001). 



10 



16 S. Izvekov and G. A. Voth, J. Chem. Phys. 116, 10372 
(2002). 

17 I-Feng W. Kuo and C. J. Mundy, Science 303 658 (2004). 

18 R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985). 

19 K. Laasonen, A. Pasquarello, R. Car, C. Lee, and D. Van- 
derbilt, Phys. Rev. B 47, 10142 (1993) 

20 J. Grossman, E. Schwegler, E. Draeger, F. Gygi and G. 
Galli, J. Chem. Phys. 120, 300 (2004). 

21 D. Asthagiri, L. R. Pratt, and J. D. Kress, Phys. Rev. E 
68, 041505 (2003). 

22 P. H.-L.Sit, M. Cococcioni, and N. Marzari, in preparation 

23 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. 
Lett. 77, 3865 (1996). 

24 N. Troullier and J. Martins, Phys. Rev. B 43, 1993 (1991). 

25 D. Vanderbilt, Phys. Rev. B 41, 7892 (1990). 

26 S. Baroni, A. Dal Corso, S. de Gironcoli, P. Giannozzi, 
C. Cavazzoni, G. Ballabio, S. Scandolo, G. Chiarotti, P. 
Focher, A. Pasquarello, K. Laasonen, A. Trave, R. Car, N. 
Marzari, A. Kokalj, http://www.pwscf.org/ 

27 D. Marx and J. Hutter "Ab-initio Molecular Dy- 
namics: Theory and Implementation", Modern 
Methods and Algorithms in Quantum Chemistry 
Forschungzentrum Juelich, NIC Series, vol. 1, (2000). 
http://www.fz-juelich.de/nic-series/Volumel/marx.pdf 

28 E. Schwegler, J. Grossman, F. Gygi and G. Galli, J. Chem. 
Phys. 121, 5400 (2004). 

29 S. Baroni, P. Giannozzi, and A. Testa, Phys. Rev. Lett. 
59, 2662-2665 (1987). 

30 M.Fuchs and M. Scheffler, Comput. Phys. 
Commun. 119, 67 (1999). Web si te: 
http : //www . f hi-berlin . mpg . de/th/f hi98md/f hi98PP 

31 Web site: http://www.nest.sns.it/~giannozz 

32 H.pbe-rrkjus.UPF and O.pbe-rrkjus.UPF 
from the vl.3 pseudopotential table at 
http : //www.pwscf . org/pseudo .htm 

33 M. Sprik, J. Hutter and M. Parrinello, J. Chem. Phys. 105, 
1142 (1996). 

34 K. Kuchitsu, Y. Morino, Bull. Chem. Soc. Jpn. 38, 805 
(1965). 

35 T. R. Dyke, K. M. Mack and J. S. Muentner, J. Chem. 
Phys. 66, 498 (1997). 

36 L. A. Curtiss, D. L. Frurip, and M. J. Blander, J. Chem. 



Phys. 71, 2703 (1979). 

37 R. M. Bentwood, A. J. Barnes, and W. J. Orville Thomas, 
J. Mol. Spectrosc. 84, 391 (1980). 

38 K. Kuchitsu, Y. Morino, Bull. Chem. Soc. Jpn. 38, 805 
(1965). 

39 R. M. Bentwood, A. J. Barnes, and W. J. Orville Thomas, 
J. Mol. Spectrosc. 84, 391 (1980). 

40 D. J. Wilber, T. Defries, and J. Jonas, J. Chem. Phys. 
65, 1783 (1976). 

41 P. Tangney and S. Scandolo, J. Chem. Phys. 116, 14 (2002) 

42 The use of norm-conserving pseudopotentials requires a 
larger plane- wave energy cutoff and thus a larger basis set. 
A smaller fictitious mass and a small time step are thus 
required. For extensive discussion, see the above reference. 

43 W. S. Price, H. Ide, Y. Arata and O. Soderman, J. Phys. 
Chem. B 104, 5874 (2000). 

44 E. Schwegler, G. Galli, and F. Gygi Phys. Rev. Lett. 84, 
2429 (2000) 

45 M. V. Ferndndez-Serra and E. Artacho cond-mat/0407237 

46 A. K. Soper, F. Bruni, and M. A. Ricci, J. Chem. Phys. 
106, 247 (1997). 

47 H. Berendsen, J. Postma, W. van Gunsteren, J. Hermans, 
In Intermolecular Forces; Pullmann, B., Ed.; Reidel: Dor- 
drecht, 1981; p 331. 

48 J. VandeVondele, F. Mohamed, M. Krack, J. Hutter, M. 
Sprik and M. Parrinello, J. Chem. Phys. 122, 014515 
(2005). 

49 D. Marx and M. Parrinello Z. Phys. B (Rapid Note) 95, 
143 (1994). 

50 D. Marx and M. Parrinello, J. Chem. Phys. 104, 4077 
(1996). 

51 M. E. Tuckerman, D. Marx, M. L. Klein and M. Parrinello, 
J. Chem. Phys. 104, 5579 (1996). 

52 D. Marx, M. E. Tuckerman, J. Hutter and M. Parrinello, 
Nature 397, 601 (1999). 

53 M. E. Tuckerman, D. Marx and M. Parrinello, Nature 417, 
925 (2002). 

54 M. Allesch, E. Schwegler, F. Gygi and G. Galli, J. Chem. 
Phys. 120, 5192 (2004). 

55 R. A. Kuharski and P. J. Rossky, J. Chem. Phys. 82, 5164 
(1985) 



